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Abstract 


Many recently developed models in areas like viscoelasticity, electrochemistry, diffu- 
sion processes, etc. are formulated in terms of derivatives (and integrals) of fractional 
(non-integer) order. In this paper we present a collection of numerical algorithms 
for the solution of the various problems arising in this context. We believe that this 
will give the engineer the necessary tools required to work with fractional models 
in an efficient way. 
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1 Introduction 


In recent years it has turned out that many phenomena in engineering, physics, 
chemistry, and other sciences can be described very successfully by models us- 
ing mathematical tools from fractional calculus, i.e. the theory of derivatives 
and integrals of fractional (non-integer) order. Some of the most prominent 
examples are given in a book by Oldham and Spanier [41] (diffusion pro- 
cesses) and the classic papers of Bagley and Torvik [2], Caputo [4], and Caputo 
and Mainardi [5] (these three papers dealing with the modeling of viscoelas- 
tic materials) as well as in the publications of Marks and Hall [34] (signal 
processing) and Olmstead and Handelsman [42] (also dealing with diffusion 
problems). More recent results are described, e.g., in the works of Chern [6] 
(finite element implementation of viscoelastic model), Diethelm and Freed 
[11] (viscoplastic material model), Gaul, Klein, and Kempfle [19] (descrip- 
tion of mechanical systems subject to damping), Glockle and Nonnenmacher 
[20] (relaxation and reaction kinetics of polymers), Gorenflo and Rutman [24] 
(so-called ultraslow processes), Metzler et al. [36] (relaxation in filled polymer 
networks), Unser and Blu [52,53] (splines and wavelets), Podlubny [44] (control 
theory), and Podlubny et al. [47] (heat propagation). Surveys with collections 
of applications can also be found in Gorenflo and Mainardi [23], Mainardi [32], 
Matignon and Montseny [35], Nonnenmacher and Metzler [39], and Podlubny 
[45]. Therefore, it seems to be fair to say that there is a significant demand 
for readily usable tools that numerically handle these mathematical problems. 
However, to date, such algorithms have been developed only to a rather lim- 
ited extent. Thus, it is the aim of this paper to identify the most important 
problems for which specific algorithms are required, to present those schemes 
that have been developed so far in a way that is directly accessible to the 
applied scientist, and to fill the remaining gaps with suitable new methods. 


The structure of this paper is as follows. In the next three sections we briefly 
recall the mathematical foundations of the fractional calculus. We then turn 
to a more detailed description and assessment of our numerical methods for 
problems described by fractional-order derivatives, integrals, and differential 
equations. After this we provide an algorithm for calculating the Mittag-Leffler 
function, which appears in solutions to fractional-order differential equations. 
Finally, after a discussion of our algorithms, we provide appendices where 
additional useful results are laid out that are connected to the numerical issues 
of fractional calculus indirectly (such as, e.g., a table of certain fractional 
derivatives, and some helpful auxiliary subroutines that can be utilized when 
working with the fundamental algorithms described in the main part of the 
paper). 
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2 Riemann-Liouville Fractional Integral 


In the classical calculus of Newton and Leibniz, Cauchy reduced the calculation 
of an n-fold integration of the function y(x ) into a single convolution integral 
possessing an Abel (power law) kernel, 

rx rx n - 1 /*xi 

J n y(x) := / / • • • / y(a? 0 ) d® 0 • • ■ dx n _ 2 da;„_ : 

Jo Jo Jo 

nsN ' i€Rt ’ (1) 

where J n is the n-fold integral operator with J° y(x) = y(x), N is the set of 
positive integers, and R+ is the set of positive reals. Liouville and Riemann 2 
analytically continued Cauchy’s result, replacing the discrete factorial (n — 1)! 
with Euler’s continuous gamma function T(n), noting that (n — 1)! = T(n), 
thereby producing [30, ^[5, Eqn. A] 

J ° y(x):= W)l \x-l'Y-° y(x ' )ix '’ a ' l€Rt ' (2) 

where J Q is the Riemann-Liouville integral operator of order a,* which com- 
mutes (i.e. J a jPy(x) = jPj a y(x) = J a+ ^ y(x) V a, 0 G R+). Equation (2) is 
the cornerstone of the fractional calculus, although it may vary in its assign- 
ment of limits of integration. In this document we take the lower limit to be 
zero and the upper limit to be some positive finite real. Actually, a can be 
complex [51], but in the applications that are of interest to us a is real. 

A brief history of the development of fractional calculus can be found in Ross 
[50] and Miller and Ross [37, Chp. 1]. A survey of many emerging applications 
of the fractional calculus in areas of science and engineering can be found in 
the recent text by Podlubny [45, Chp. 10]. 


3 Caputo-Type Fractional Derivative 

3.1 Fundamental definition 


From this single definition for fractional integration, one can construct several 
definitions for fractional differentiation (cf., e.g., with Refs. [45,51]). The spe- 
cial operator that we choose to use, which requires the dependent variable 

2 Riemann’s pioneering work in the field of fractional calculus was done during 
his student years, but published posthumous — forty-four years after Liouville first 
published in the field [50]. 
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y to be continuous and [a] -times differentiable in the independent variable x, 
is defined by 

D° y{x) := jN-^N y ( x ), (3) 

such that 

Jim D“ y(x) = D n y(x) for any n e N, (4) 

with D®y(x) = y(x), where [a] is the ceiling function giving the smallest 
integer greater than or equal to a, and where a — >■ n_ means a approaches n 
from below. The operator D n , ns N, is the classical differential operator. It is 
accepted practice to call £>“ the Caputo differential operator of order a, after 
Caputo [4] who was the among the first to use this operator in applications 
and to study some of its properties. 3 Appendix A presents a table of Caputo 
derivatives for some of the more common mathematical functions. 

The Caputo differential operator is a linear operator, i.e. 

D * (ay + fez) (x) = a£>? y(x) + fe£>* z(x) (5) 

for arbitrary constants a and fc, that commutes, viz. 

D;v?y(x) = D f .D?y(x) = IK +e y(x) Va,0eR+, (6) 

if y(x) is sufficiently smooth, and it possesses the desirable property that 

D°c = 0 for any constant c. (7) 

The more common Riemann-Liouville fractional derivative D a , although lin- 
ear, need not commute [45, pg. 74]; furthermore, D a c = dHjW'“c = 
crr _Q /r(l — a), which depends on x. Ross [50] attributes this startling fact 
as the main reason why the fractional calculus has historically had a difficult 
time being embraced by the mathematics and physics communities. 


3 Actually, Liouville introduced the operator in his historic first paper on the topic 
[30, f6, Eqn. B], Still, nothing in Liouville’s works suggests that he ever saw any 
difference between D° := jbd -0 !)! 0 ! and D a ■= pW jf Q l _a ) with D a being his 
accepted definition [30, first formula on pg. 10] — the Riemann-Liouville differential 
operator of order a. Liouville freely interchanged the order of integration and differ- 
entiation, because the class of problems that he was interested in happened to be a 
class where such an interchange is legal, and he made only a few terse remarks about 
the general requirements on the class of functions for which his fractional calculus 
works [31]. The accepted naming of the operator D° after Caputo therefore seems 
warrented. 

Rabotnov [49, pg. 129] introduced this same differential operator into the Russian 
viscoelastic literature a year before Caputo’s paper was published. Regardless of this 
fact, operator D ® is commonly named after Caputo in the present-day literature. 
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The Riemann-Liouville integral operator J a and the Caputo differential oper- 
ator D“ are inverse operators in the sense that 


D+J a y{x) = y(x) and J a D° y{x) = y(x) - £ y$, a£l+, (8) 

with yj,+ := D k y(0 + ), where |_aj is the floor function giving the largest integer 
less than or equal to a. The classic n-fold integral and differential operators 
of integer order satisfy like formulae, viz.: D n J n y(x) = y(:r), and J n D n y(x) -- 
y(x) — o yj+, n € N. In this sense, the Caputo derivative and Riemann- 
Liouville integral are analytic continuations of the well-known n-fold derivative 
and integral from the classical calculus. 

Remark 1 Fractional derivatives do not satisfy the Leibniz product rule as it 
is known in integer-order calculus. For example, whenever the Caputo deriva- 
tive is restricted so that 0 < a < 1, its Leibniz product rule is given by 

n?(y * ■)<*> = rfS * + (^)<*> * ■<*> 

+ £ (y (^“y)M x (D*z) (re), (9) 

wherein, unlike the Leibniz product rule for integer- order derivatives, the bi- 
nomial coefficients (fy = a(a — l)(o: — 2) • • • (a — k + l)/fc! (with = 1, 
a £ R+, and k e N) do not become zero whenever k > a because a £ N 
(i.e., the binomial sum is now of infinite extent). A similar infinite sum exists 
for the Leibniz product rule of the Riemann-Liouville fractional derivative ( cf. 
Podlubny [45, pp. 91-97]). 


3.2 Alternative integral expressions 


The Caputo derivative defined in Eq. (3) can be expressed in a more explicit 
notation as the integral 


***«> = f ( , - ^-w ( oro Mw d *'- do) 


where the weak singularity caused by the Abel kernel of the integral operator 
is readily observed. This singularity can be removed through an integration 
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by parts 


D a v( x) — ( rM-° V C Q 1) 

* yw r(i + h - a) \ y ° + 

+ (x - X ’) M “(D 1+ r«ly)(i')dx'V (11) 

provided that the dependent variable y is continuous and (1+ [a])-times dif- 
ferentiable in the independent variable x over the interval of differentiation 
(integration) [0, a?]. In Eq. (11) the power-law kernel is bounded over the en- 
tire interval of integration; whereas, in Eq. (10) the kernel is singular at the 
upper limit of integration. 

The two representations in formulas (10) and (11) are quite useful for pen- 
and-paper calculations, but in order to obtain a numerical scheme for the 
approximation of such fractional derivatives, we found it even more helpful to 
look at yet another representation that seems to have been introduced into 
this context by Elliott [14]; namely, 

D ° y(x) = rfby i: (X - y(x ' ) d1 '’ < 12 > 

This representation can also be obtained from Eq. (10) using the method 
of integration by parts, but with the roles of the two factors interchanged. 
The advantage here is that the function y itself appears in the integrand in- 
stead of its derivative. The disadvantage is that the singularity of the kernel 
is now strong rather than weak, and thus we have to interpret this integral as 
a Hadamard-type finite-part integral. This is cumbersome in pen-and-paper 
calculations but, as we shall see below, it is not a problem to devise an al- 
gorithm that makes the computer do this job. We provide a brief description 
of such an algorithm in the following pages. For more details, the interested 
reader is referred to Refs. [7,14] and the references cited therein. 


4 Caputo-Type FDE’s 


Fractional material models, a topic of interest to one of the authors (ADF), 
are an important application where systems of fractional-order differential 
equations (FDE’s) arise that need to be solved in accordance with appropriate 
initial and boundary conditions. A number of other examples are well known; 
we refer the reader to the book by Podlubny [45, Chap. 10] and the survey of 
Mainardi [32] for more information. 
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A FDE of the Caputo type has the form 

D°y{x) = f (x, y(x)) , a, x G !+, (13) 

satisfying a set of (possibly inhomogeneous) initial conditions 

D fc y(o + ) = yj+. k = 0,l,...,[a\, (14) 

whose solution is sought over an interval [0, X], say, where X G R+. It turns 
out that under some very weak conditions placed on the function f of the 
right-hand side, a unique solution to Eqs. (13 & 14) does exist [8]. 

A typical feature of differential equations (both classical and fractional) is the 
need to specify additional conditions in order to produce a unique solution. 
For the case of Caputo FDE’s, these additional conditions are just the static 
initial conditions listed in (14), which are akin to those of classical ODE’s, 
and are therefore familiar to us. In contrast, for Riemann-Liouville FDE’s, 
these additional conditions constitute certain fractional derivatives (and/or 
integrals) of the unknown solution at the initial point x = 0 [27], which are 
functions of x. These initial conditions are not physical; furthermore, it is not 
clear how such quantities are to be measured from experiment, say, so that 
they can be appropriately assigned in an analysis. 4 If for no other reason, the 
need to solve FDE’s is justification enough for choosing Caputo’s definition 
(i.e. , D“ := for fractional differentiation over the more commonly 

used (at least in mathematical analysis) definition of Liouville and Riemann 
(viz., D a := HWjM-“). 


5 Numerical Approximation of Caputo-type Derivatives 

5.1 Fundamental ideas 


Unlike ordinary derivatives, which are point functionals, fractional deriva- 
tives are hereditary functionals possessing a total memory of past states. A 

4 We explicitly note, however, the very recent paper by Podlubny [46] who attempts 
to give highly interesting geometrical and physical interpretations for fractional 
derivatives of both the Riemann-Liouville and Caputo types. These interpretations 
are deeply related to the questions: What precisely is time? Is it absolute or not? 
And can it be measured correctly and accurately, and if so, how? Thus, we are 
still a long way from a full understanding of the geometric and physical nature 
of a fractional derivative, let alone from an idea of how we can measure it in an 
experiment, but our mental image of what fractional derivatives and integrals ‘look 
like’ continues to improve. 
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numerical algorithm for computing Caputo derivatives has been derived by 
Diethelm [7] 5 and is listed in Alg. 1. Validity of its Richardson extrapolation 
scheme for 1 < a < 2, or one similar to it, has to date not been proven, 
or disproven. Here y n denotes y(x n ), while y N represents y(X) where [0, X] 
is the interval of integration (fractional differentiation) with 0 < x n < X. 
This algorithm was arrived at by approximating the integral in Eq. (12) with 
a product trapezoidal method, thereby restricting 0 < a < 2. Similar algo- 
rithms applicable to larger ranges of a can be constructed by using the general 
procedure derived in Ref. [7], if they become needed. 


Algorithm 1 

Computation of a Caputo fractional derivative (0 < a < 2, a ^ 1): 

For interval [0, X ] with grid {x n = nh : n = 0, 1, 2, ... , N} where h = X/N , 
compute 


D?y N (h) = [ft-“/r(2 - a)] E" =0 «„,jv (y«-„ - e£ J „[(W - n) k h k /k\ Jyj?) 

o;y(x) = o;y J v(A) + 0(ft 2 -), 

using the quadrature weights (derived from a product trapezoidal rule) 

1, if n — 0, 

&n,N — < (n + l) 1-a — 2n 1-Q + (n — l) 1-a , if 0 < n < N, 
^(l-a)N~ a - N 1 - 0 + (N -l) 1 - 0 , if n = N. 

Refine, if desired, using Richardson extrapolation 


o; y ;(x) = (D; y ;:l(x) - r— o; y r‘(^))/( i - 2 "—), 
o; y (x) = D; y ;(x) + o(h'-), 

such that if 0 < a < 1 then r u _i is assigned as 


> 


r 0 = 2 — a, 
n = 2, r 2 = 3 - a, r 3 = 4 - a, 
rA = 4, r 5 = 5 - cc, r 6 = 6 - a, 
r 7 = 6, • • • . 


The Griinwald-Letnikov algorithm is often used to numerically approximate 
the Riemann-Liouville fractional derivative (cf., e.g., with Oldham and Spanier 
[41, §8.2] and Podlubny [45, Chp. 7]), and it was the first algorithm to appear 
for approximating fractional derivatives (and integrals). 


5 Apparently this algorithm first appeared in the PhD thesis of Chern [6], unbe- 
knownst to us (KD) at the time of writing Ref. [7]. Chern used this algorithm to 
differentiate a Kelvin- Voigt, fractional-order, viscoelastic, material model in a finite 
element code. He did not address stability or uniqueness of solution issues; he did 
not compute error estimates; and he did not utilize an extrapolation scheme to 
enhance solution accuracy. 
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Diethelm’s Quadrature Weights for Fractional Differentiation 

x/X, 0<x<X 



Fig. 1. Weights of quadrature a n> jv for approximating Caputo’s fractional derivative 
(12) over interval [0, X] using Diethelm’s [7] Alg. 1, plotted here for various values 
of a with N = 8. 


5.2 Quadratures 


The extent of rememberance of past states exhibited by the hereditary nature 
of a fractional derivative is manifest, for example, in its weights of quadrature, 
as illustrated in Fig. 1. This operator exhibits a fading memory: 0.001 < 

1 08,8 1 < 0.01 for the six cases plotted in this figure. If Dy(X) were to be 
approximated by a backward difference with h = X/8, then the effective 
weights of quadrature would be Oo,8 = 1 and ai >8 = — 1 with all remaining 
weights being zero, as represented by the line segments in this figure. Similarly, 
if D 2 y(X) were to be approximated by a like backward-difference scheme, then 
a 0 ,8 = 1, a 1>8 = —2 and 02 )8 = 1 with all remaining weights being zero. It is 
evident from the data presented in Fig. 1 that the weights of quadrature a„ )8 
for approximating D°y(X) are compatible with those for the first- and second- 
order backward differences, and that fractional quadratures have additional 
contributions that monotonically diminish with increasing nodal number from 
node n = 2 fading all the way back to the origin at node n = N. This suggests 
that a truncation scheme may be able to be used to enhance algorithmic 
efficiency for some classes of functions, but not all. 
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5.3 Richardson extrapolation 


Richardson extrapolation is a technique that can often be used to increase the 
accuracy of results [12], As we utilize it, this technique follows a triangular 
scheme — a Romberg tableau — that has the form 

d; y 0 ° 

d; y? d; yi 

D; y\ DZyl D? y\ (15) 

D!y °3 Dty\ d ; yf c; y | 


Constructing the first column of the tableau constitutes the bulk of the com- 
putational effort. Here D° y° := -D“yw(/i) denotes the value of D “ y evalu- 
ated numerically at X over [0, AC] using an initial stepsize of h (= X/N ), 
D+y® '■= D+y N {h/2) is computed using the refined stepsize of \h (= X/2N), 
Dy° 2 := D*y N (h/4) is computed using the even more refined stepsize of \h 
(= X/4N), while y° := D°y^(h/ 8) is computed using the further refined 
stepsize of | h (= X/8N), etc. The remaining columns are quickly evaluated 
using the recursive formula listed in Alg. 1. The advantage of constructing this 
tableau is that the elements D“ y“(X) in the w th column converge for fixed u 
and increasing v towards the true value of the Caputo derivative as 0(h ru ). 
Hence, the further one moves to the right in the tableau the faster the column 
converges, and this level of convergence requires less computational effort to 
achieve than a direct computation of D°y^(X-, h ) when computed to a similar 
accuracy of O(h ro ) ~ 0(h ru ). 


5.4 Step-size choice 


The error analysis mentioned above is only a truncation error analysis. It 
assumes that the calculations are done in exact arithmetic, and it does not 
take into account effects like roundoff. When one needs to look at these effects 
too, it is possible to ask for a step size h = X/N whose combined effect arising 
from both error sources is minimized. As we have seen above, it is likely 
that the truncation error decreases with the step size h, whereas roundoff 
tends to have the opposite behavior, so we should be looking for a sort of 
compromise. The considerations in this context are very similar to those for 
integer-order derivatives [48, §5.7]. Roughly speaking, it turns out that the 
roundoff error behaves as h~ a e y y(ri), where e y is the relative accuracy with 
which one can compute y, and where r\ designates some number within the 
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y = o“x 



Fig. 2. The Caputo derivative D%x for various values of a. 


interval [0, X]. Moreover, the truncation error is close to c Q /i 2_Q! y" (£) , where c a 
is some constant independent of y, and £ is some other number also contained 
in the interval [0, A - ]. Consequently, an optimum step size would be of order 
h ~ fa/y^/y^f)] 1 / 2 when minimizing with respect to both trucation and 
roundoff errors. 

Unless specific information indicating the contrary is available, one may as- 
sume that y and y" are not too irregularly behaved. Under these conditions 
y(v) ~ y(X) and y"(r)) ~ y"(X), and one can then follow the suggestion 
of Press et al. [48, p. 187] by setting y(X)/y"(X) « X (except near X = 0 
where some other estimate for this quantity should be used). This scheme pro- 
vides some advice on the choice of step size if roundoff effects are considered 
problematic in some specific application at hand. 


5. 5 Examples 


The simplest fractional-order derivative that one can consider is D“y(x) where 
y(x) = x. A plot of this derivative is presented in Fig. 2. There are three special 
cases shown in this figure: D°x = x, D\x = 1, and D*x = 0 V a > 1, which are 
in accordance with the classic definition for differentiation. The intermediate 
values of a = 1/4, 1/2, and 3/4 suggest a smooth transition in the response 
between a = 0 and 1, as one would expect. 
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0 


- 0.2 


- 0.4 


- 0.6 


- 0.8 


x 

Fig. 3. The Caputo derivative D" e~ x for various values of a. 

The second example, presented in Fig. 3, is for D“y(x) where y(x) — e~ x . 
There is a discontinuous jump at x = 0 between the first-order derivative 
De~ x = — e~ x and its Caputo counterparts D“e -X given that 0 < a < 1. 
These Caputo derivatives start at zero, drop rapidly piercing the — e - * curve, 
and then gently swing up eventually asymptoting to zero from below. 


6 Numerical Approximation of Riemann-Liouville Integrals 

6.1 Fundamental ideas 

Over the course of our research we have not only had to approximate frac- 
tional derivatives, but also fractional integrals. As indicated above, the natu- 
ral concept for the fractional integral to be used in connection with Caputo 
derivatives is the Riemann-Liouville integral described in Eq. (2). We there- 
fore present a numerical scheme for the solution of this problem, too. The 
underlying idea of the algorithm, stated in a formal way in Alg. 2 below, is 
completely identical to the idea presented above for the Caputo derivative; 
that is, we employ a product integration technique based on the trapezoidal 
quadrature rule. Said differently, we replace the given function y(x') on the 
right-hand side of Eq. (2) by a piecewise linear interpolant, and then we cal- 
culate the resulting integral exactly. As a matter of fact, this algorithm will 
also be part of the scheme introduced in Alg. 3 in the pages that follow for the 
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numerical solution of a certain type of fractional-order differential equation. 

Algorithm 2 

Computation of a Riemann-Liouville fractional integral (a > 0): 

For interval [0, X] with grid {x n = nh : n — 0, 1, 2, , N} where h = X/N, 
compute 

J a y N (h) = [h a /r(2 + a)] EiJLo Cn,N y n , 

J a y(X) = J a y N (h) + 0(h 2 ) , 

using the quadrature weights ( derived from a product trapezoidal rule ) 

\1 + a)N a - N 1+a + (N- l) 1+a , if n = 0, 

c n ,N = <(A-n+l) 1+a - 2 (N-n) l+a + (N-n- 1) 1+Q , if 0 < n < N, 

1, if n = N. 

Refine, if desired, using Richardson extrapolation 

J° yi(X) = (j°ytl(X) - 2 J''yr'<X>)/n - 2 '—) , 

J«y(AT) = J“y-(X) + 0(h'~), 
such that if0<a<l then r u _i is assigned as 

r 0 = 2, r x — 2 + a, 

= 3, r 3 = 3 + a, 

r 4 = 4, r 5 = 4 + a, 

r 6 = 5, • • • . 

Note 1 Whenever a > 1 in Alg. 2, the same values appear in the sequence 
r 0 i r ii r 2 i ■ ■ ■> but they must be re-ordered to keep the sequence monotonic. For 
example, if 1 < a < 2 then we have r 0 = 2, rq = 3, r 2 = 2 + a, r 3 = 4, r 4 = 
3 + a, — 

It is easily seen that the error of this algorithm is of the order 0(h 2 ) where, 
as in Alg. 1, h denotes the step size. Once again, we can improve the accuracy 
by adding a Richardson extrapolation procedure to the plain algorithm. The 
required exponents are known (cf. Ref. [26, §4]) and the resulting scheme is 
detailed in Alg. 2. Both the fundamental algorithm itself, and the Richardson 
extrapolation procedure, may be used for any positive value of a; there is no 
need to impose an upper bound on the legal range for a. This is due to the 
fact that the Abel (power law) kernel in the definition (2) of the Riemann- 
Liouville integral is regular, or at worst, weakly singular, and hence, integrable 
(at least in the improper sense) for any a > 0. In contrast, the corresponding 
kernel in representation (12) of the Caputo derivative is not integrable. This 
kernel requires special regularization methods that are compatible with our 
approximation method, and as such, our scheme for approximating Caputo 
derivatives is only valid for 0 < a < 2; whereas, our corresponding scheme for 
approximating Riemann-Liouville integrals is valid for all a > 0. 
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Table 1 


Averaging procedure used to compute effective weights of quadrature for approxi- 
m ating Riemann-Liouville integration as they relate to Alg. 2. 


Subinterval 

Averaged Quadrature Weight (c n , jv) 

[0,X/N] 

(co,n) = Co, AT + \ Cl, AT 

[nX/N, (n + 1)X/N] 

( c n,iv) = 2 { c n>N + c n+l,iv) ? (?T = 1 , 2, . . . } N — 2) 

[(N-1)X/N,X) 

( CN-1,N ) = \ Cat-1, AT + CN,N 


Notice the formal correspondence between Alg. 2 (for fractional integration of 
order a) and Alg. 1 (for fractional differentiation of order a). Except for the 
initial conditions that have to be taken into account additionally, the latter is 
simply obtained from the former by replacing the parameter a by —a. 6 This 
relates to the intuitive (but not mathematically strictly correct) interpretation 
of fractional differentiation and integration being inverse operations. Also no- 
tice that the index ordering is inverted between these two algorithms, which 
is in keeping with accepted indexing conventions. Algorithm 1 indexes from 
x 0 = X to xn = 0, while Alg. 2 indexes from x 0 = 0 to = X. 


6.2 Quadratures 


A visualization of quadrature weight versus nodal index for several values of 
a pertaining to Alg. 2 is presented in Fig. 4. What is striking about this 
figure is the obvious difference between domains 0 < a < 1 and 1 < a < 2. 
Whenever a = 1, the algorithm reduces to classic trapezoidal integration. 
Whenever 0 < a < 1, the earlier states will contribute less to the overall 
solution than will the more recent states, but they do not entirely fade out. 
Fractional integration exhibits a long-term memory loss when 0 < a < 1 but, 
unlike fractional differentiation, fractional integration does not experience a 
rapid loss (or fading away) of past memories. Also, the smaller the value of a 
(i.e., the closer it is to zero) the greater the degree of long-term memory loss 
will be. In contrast, whenever 1 < a < 2, the earlier states will contribute more 
to the overall solution than will the more recent states. Fractional integration 
therefore exhibits a short-term memory loss when 1 < a < 2. This is like an 
elderly person who remembers in vivid detail what happened years ago, but 
who cannot recall what took place yesterday. Furthermore, the greater the 
value of a (i.e., the closer it is to two) the more pronounced the short-term 
memory loss will be. 


6 Similarly, the Griinwald-Letnikov algorithm for approximating Riemann-Liouville 
fractional derivatives of order a also applies for approximating Riemann-Liouville 
fractional integrals by replacing their algorithmic parameter a with -a [41, §8.2], 
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Diethelm’s Quadrature Weights for Fractional Integration 

x/X, o < x < x 

0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1 



Fig. 4. Effective weights of quadrature (c n) jv) for approximating the Rie- 
mann-Liouville fractional integral (2) over interval [0, X] using Alg. 2, plotted here 
for various values of a with N — 8. 

The line segments displayed in Fig. 4 represent normalized weights of quadra- 
ture averaged over each subinterval. The actual nodal weights, c ni n, are often 
observed to be non-monotonic at either of the two nodal endpoints. In this 
integration scheme there are N + 1 nodal weights that apply to N subin- 
tervals, but there should be exactly one weight per subinterval. So how the 
algorithm works (internally, and roughly speaking) is to average these weights 
in a trapezoidal fashion, as outlined in Table 6.2. In other words, the inner 
weights are divided into two equal halves with each half going to one of the 
two adjoining subintervals. In addition to averaging, the displayed line seg- 
ments have been normalized over the interval [0, 1]. Normalizing allows one to 
discern the influence of a on quadrature in a meaningful way. The outcome 
is an averaged and normalized quadrature weighting that is monotonic in the 
nodal index number, as demonstrated by the line segments in Fig. 4, where 
there is a monotonic increase (decrease) in the effective weight of quadrature, 
( c n,Af) j ma x m (c m)JV ), with increasing nodal index number for 0 < a < 1 
(1 < a < 2). 

4 

6.3 Optimizing performance 


In a recent paper, Ford and Simpson [16] developed a faster scheme for the 
calculation of fractional integrals. A reduction in the amount of computational 
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work can be achieved by using a graded mesh, thereby taking the 0(N 2 ) 
method of Alg. 2 and making it 0(N log N). The underlying idea is based on 
the fact that the fractional integral possesses a fundamental scaling property 
that can be exploited in a natural way. 


If q > 0 and the function y is suitably integrable then, with J a defined as 
before, we observe that 


J°y(qx) 



yW 

(qx — x') l ~ a 


dx'. 


We make the substitution of x' by qx' to obtain 


J a y(qx) 


<f [ x y(ffc') i / 
r(a)J 0 (x-x'y-° ’ 


to which it follows that, for any p G N, 


J a y(q p x) 


qpa I>x 

r(a) Jo 


ytfV) 

(x — x') 1 ~ a 


dx'. 


(16) 


(17) 


(18) 


Thus, the weights necessary to calculate Cl%y(nh) « J a y(nh) (using a step 
length of h) can be used to calculate £lq Ph y(nq p h) « J a y(nq p h) (using a step 
length of q p h) for any p G N simply by multiplying the resulting sum by q pa . 
That is, 

fifty ( nh ) = 12 u n-jy(jh) O fiq Ph y(nq p h) = q pa ^n-jy(jq p h). (19) 

j = o 3=0 

We use a sequence of nested meshes in the following way: for h G M+, we 
define the mesh M/, by Mh — { hn , n G N}. If p, q, r G N, q > 0, r > p, then 
MqPh D M q r h . We then decompose the interval [0,x], for fixed X > 0, in the 
following way: 

[0, x] = {0,x-q m X}u[x-q m X,x-q m - 1 X}U---U[x-qX,x-X}u[x-X ) x}, 

where m G N is the smallest integer such that x < q m+1 X. 

Now, as is clear from the definition of J a , the singularity in the kernel of 
the fractional derivative (integral) occurs when x = x', while the value of the 
kernel l/(x - x') l ~ a ->■ 0 as x oo for x' -4 0. We combine this insight with 
the scaling property, and we use a step length of h over the most recent time 
interval [x - X, x\ and successively larger step lengths over earlier intervals 
in the following way: let x, X, h G R, q m+1 X > x > q m X, X > 1, h > 0 
with x — nh for some n G N. We can therefore rewrite the Riemann-Liouville 
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Fig. 5. The Riemann-Liouville integral J°x for various values of a. 
fractional integral in a form suitable for efficient computation; it being, 

m — 1 

J [My(*) = J{ X -x,x\y( x ) + E + ^o,x- 9 -x]y(^) 

2 — 0 

771—1 

= J [x-X,x\y(x) + E 9 <QJ [x- 9 *+^,x- 9 *x]y(9^) + q ma J^ gmx] y(q m x). (20) 

i=0 

This reduces the computational effort required, both to calculate the weights 
we use (since now we need fewer different weights) and (more importantly) to 
evaluate the integral. We employ the same weights that are listed in Alg. 2, 
but the above scheme is independent of the algorithm used for integration. 


6.4 Examples 


To complement Fig. 2, consider J a y(x) where y(x) = x, which is plotted in 
Fig. 5 for values of a that vary between 0 and 2. The well-known solution 
J x x = |x 2 bisects the displayed family of curves, while the fractional (in a) 
solutions smoothly traverse between the integer solutions. 

Our second example for fractional integration is akin to our second example 
for Caputo differentiation, viz., J a y(x ) where y(x) = e -x , which is plotted in 
Fig. 6. The curves in this figure are essentially mirror images of those presented 
earlier in Fig. 3 for D“e -X , except that the a-ordering is reversed. 
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x a -x 

J e 



Fig. 6. The Riemann-Liouville integral J a e x for various values of a. 

7 Numerical Approximation of Caputo-Type FDE’s 


7.1 A basic general-purpose algorithm 


A numerical algorithm that solves Caputo-type FDE’s has been derived by 
Diethelm et al. [10] and is listed in Alg. 3. A thorough analysis of its algo- 
rithmic error is given in [9] . This algorithm is of the PECE (Predict-Evaluate- 
Correct-Evaluate) type. Other numerical algorithms exist that solve FDE’s 
(e.g., Gorenflo [21] and Podlubny [45, Chp. 8]), but they focus on solving 
Riemann-Liouville FDE’s and usually restrict the class of FDE’s to be linear 
with homogeneous initial conditions. Algorithm 3 solves non-linear Caputo 
FDE’s with inhomogeneous initial conditions, if required. 

The restriction that a 7^ 1 in this algorithm is purely for formal reasons. If 
a = 1, then we can still implement the algorithm exactly in the indicated 
way. It must be noted, however, that it then is the limit case of an algorithm 
for fractional differential equations, and these equations involve non-local dif- 
ferential operators. Thus, the resulting scheme is non-local, too. In contrast, 
a method constructed specifically for a first-order equation will, in practice, 
always make explicit use of the local structure of such an equation to save 
memory and computing time. Therefore, the case a = 1 of our algorithm 
will never be a competitive alternative to the usual methods for first-order 
equations. In particular, our algorithm is distinct from the algorithm known 
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as the second-order Adams-Bashforth-Moulton technique for first-order prob- 
lems when a: = 1. 


Algorithm 3 

Computation of a Caputo FDE (0<a<2,a^l): 

For interval [0, X] with grid {x n = nh: n = 0, 1, 2, . . . , N: h = X/N}, 
predict with 

y 5(A) = + [A“/r(i + <»)] ES K» f(*», y») , 

using the quadrature weights (derived from, a product rectangular rule) 
b ; N = (N-n)°-(N-n- 1)» , 
and erjafuate f(X, yjy), then correct with 

Y N(h) = Z\?J 0 (X k /k\)y$ 

+ [h a /T(2 + a)] (T,n=o °n,N f(®„, y„) + cyv.jv f(A, y£)) , 
y{X) = y N (h) + 0(h min ( 1+a >V) , 

using the quadrature weights (derived from a product trapezoidal rule) 

'(1 + a)N a - N l+a + (N - l) 1+a , if n = 0, 

c n ,N = <(7V— n+l) 1+Q - 2(A-n) 1+Q + (N-n- 1) 1+Q , if 0 < n < N, 

U, if n = N, 

and finally, re- ermfaate f(X, y^) saving it as f(x^,y n)- 
Refine, if desired, using Richardson extrapolation 

y;w = (y":!W - 2 r ~‘yr 1 W)/U - 2 '—) , 
yW = y«W + o(/i r “) , 

such that whenever 0 < a < 1 the exponent r u _i is assigned as 

r 0 = 1 + a, 
r i = 2, r 2 = 2 + 0 ;, r 3 = 3 + a, 
r 4 = 4, r 5 = 4 -I- a, r 6 = 5 + a, 
r 7 = 6, • • • , 

or whenever 1 < a <2 it is assigned as 
r 0 = 2, ri = 1 -(- a, r 2 = 2 + a, 

^3 = 4, r 4 = 3 + a, r 5 = 4 + a, 
r 6 = 6, • • • . 


7.5 Quadratures 

Illustrations of quadrature weight versus nodal index for several values of a, 
as they pertain to the PECE method of Alg. 3, are presented in Figs. 4 & 7 


19 



Diethelm’s Quadrature Weights for Fractional Integration 
x/X, 0 < x < X 

0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1 



Fig. 7. Normalized weights of quadrature b n ^ for the predictor that approximates 
Caputo FDE’s (13) over interval [0,X] when using Diethelm et al.’s [10] Alg. 3, 
plotted here for various values of a with N = 8. 

with the former figure pertaining to the corrector and the latter pertaining to 
the predictor. FDE’s, like fractional integrals, exhibit long-term memory loss 
when 0 < a < 1, no memory loss when a = 1, and short-term memory loss 
when 1 < a < 2. 

Unlike the (c n<N ) in Fig. 4, where the N + 1 quadrature weights are averaged at 
the beginning and end of each subinterval in order to get N effective weights 
for N subintervals, the b n>N in Fig. 7 are fixed to the beginning of each of 
the N subintervals, and as such, do not require any ‘effective averaging’ to 
take place. This is a consequence of the b n , N quadrature weights belonging to 
an explicit integrator, while the c n>N weights belong to an implicit integrator. 
Contrasting Figs. 4 & 7, there is little difference between the b n< N and (c n! jv) 
curves, indicating that there is a much stronger influence of a on the weights 
of quadrature than there is on the order of accuracy (e.g., 0(/7 mm ^ 2,1+Q ^)) that 
a particular integration scheme provides. 

Quadratures for non-equidistant nodes have also been derived by Diethelm 
and Freed [11], but they are vastly more expensive to compute, and as such, 
are not the preferred quadratures to use. 

Differential equations of fractional order have found recent applications in a 
variety of fields in science and engineering (e.g., see references in [27,32,45]): 
chemical kinetics theory; electromagnetic theory; transport (diffusion) theory; 
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fractals, splines, and wavelets; control theory; electronic circuit theory; porous 
media; etc. 


7 . 3 Efficient approximations 


The nested mesh algorithm for the evaluation of fractional integrals given by 
Ford and Simpson [16] and described in §6.3 can also be adapted to provide 
an effective fast algorithm for the solution of Caputo FDE’s. Whereas, the 
unmodified method possesses O(N) operation counts at each stage and 0(N 2 ) 
overall, this modification leads to a more efficient scheme with 0(N log N) 
counts overall, while retaining the accuracy of the method. 


7.4 Examples 


Many physical processes can be represented by a quantity that begins at some 
nominal value, and later asymptotes to some saturated state under steady 
conditions, e.g., we [11] have used such an equation to describe the evolution 
of back stress (an internal state variable) in viscoplasticity. This effect can 
be modeled by the differential equation D+y(x) = 1 — y(x) whose solution is 
presented in Fig. 8 for values of a less than 1, and in Fig. 9 for values of a 
greater than 1. The solution will be monotone provided that 0 < a < 1; it 
will be a damped oscillator whenever 1 < a < 2; it will be a perfect oscillator 
if a = 2; and it will become an unstable oscillator for any value of a > 2. 


8 Mittag-Leffler Function 


Analytic solutions to fractional-order differential equations are often expressed 
in terms of the Mittag-Leffler function, and so we also need a scheme that 
computes this special function, which is the topic of this section. 


8.1 Fundamental definitions 


The (generalized) Mittag-Leffler function E a ^(z) is an entire function (in z e 
C) of order l/a that is defined by the power series [15, §18.1] 


EcAz) = £ 


ti r (/9 + ak) • 


ot £ K+ , /? £ R, z £ C, 


( 21 ) 
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Fig. 8. The solution of the Caputo differential equation D“y(x) = 
various values of a € (0, 1] satisfying an initial condition of y(0) = 0. 


D“y(x) = 1 - y 

y(0) = 0, Dy(0) = 0 

2 1 1 [ 1 I ' T 


1.5 


a= 1 

a =1.25 

a= 1.5 

a= 1.75 



1 — y(x ) for 
»y(0) = 0. 









whose derivative is 


pi r~\ — dE a> p(z) _ (1 + k)z k 

“ J,W_ dz £;r[/J + a(l + *)]’ 

with E a (x ) := E a> i(x) being the original function studied by Mittag-Leffler 
[38]. This function plays the same role in differential equations of fractional 
order that the exponential function e 2 plays in ordinary differential equations; 
in fact, Ei t i(z) = e 2 . 

A special form of the Mittag-Leffler function, 

G(t - t') := J5 ai i(— [(t - t')/r] a ^, 0 < a < 1, 0 < r, 0 <t'<t, 

and its derivative, 

M(t - t') - ~ ^1 _ -ffq,oH(*-* , )/' r ] Q ) 

{ 1 ' dt' t - V 

appear in the theory of fractional-order viscoelasticity (FOV) [17]; it is also 
of major significance in any application using linear FDE’s with constant 
coefficients [23,24]. 

We now present some important properties of the Mittag-Leffler function, and 
a numerical algorithm for its accurate evaluation, both of which are useful 
when considering differential equations of fractional order. 



8.2 Analytical properties 


In spite of the fact that in applications to differential equations of fractional 
order the Mittag-Leffler function is typically restricted to the real line, we still 
need to give some of its properties in the complex plane. The main reason 
for this is that the numerical algorithm we present in the following pages 
consists of two parts: the first part gives a numerical value for the Mittag- 
Leffler function with a < 1, while the second one uses, for a > 1, some special 
formulae that reduce this case to the previous one. These special formulae are 
defined over the complex plane and are given by 


Ea,p(z) 

and 


-| m 

^/Pm + i ) ^(^ 1/<2m+1) e‘ 2 '‘ /|2m+,) ) . rn = 0,1,2 

(23) 


h— — m 


-J 771—1 

= - E E alm A^ lm ^’ hlm ), m = 1,2 

’ll I A 


(24) 


h=0 
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where i (:= \f—i) is the imaginary unit number. Obviously, if a > 1, and even 
if 2 is a real number, we still need to evaluate the numerical values of the 
Mittag-Leffler function with, in general, a complex argument (and a < 1). 

First, we present some important integral representations of the Mittag-Leffler 
function. Let us denote by 7 (p; <p) (p > 0, 0 < <p < n) a contour in the complex 
A-plane with non-decreasing argA consisting of the following parts: 

(1) the ray argA = — <p, |A| > p; 

(2) the arc — ip < argA < <p from the circumference |A| = p; 

(3) the ray argA = (p, |A| > p. 

In the case where 0 < ip < n, the complex A-plane is divided into two un- 
bounded parts by the contour 7 (p; ip): domain <p) is to the left of the 

contour, while domain G^(p; <p) is to its right. If <p = n, the contour 7 (p; <p) 
consists of the circumference |A| = p and the cut — 00 < A < —p. In this 
case, domain G( _ )(p;<p) becomes the circle |A| < p, while domain G^ + \p\p) 
becomes the region {A : | arg A| < n, |A| > p}. 

Let 0 < a < 2, let (3 be an arbitrary (real or complex) number, and let a 
non-negative number <p be chosen such that 

an . , , , x 

— < <p < mm{7r, cc7r}. (25) 

A 

Then we have the following integral representations for the Mittag-Leffler 
function: 


Ea,p(z) 


1 r e Al/ “ At 1 -^ 
i27TO: Jy(p-,ip) A — z 


dA, zeG ( \p-,v), 


(26) 


and 

Jl-p)/a e z l / a 

E Ql p(z) = h 

a 


[ 

Vl^KOt Jy(p;<p) 


e * 1/a A (i->9)/« 
A — z 


dA. zetf+Hfip). (27) 


If f3 is a real number, as assigned in Eq. (21), then the formulae of Eqs. (26 & 
27) can be rewritten in forms that are more suitable for numerical evaluation 
(cf. Gorenflo et al. [22]). In particular, if 0 < a < 1, /3 G R, | arg z | > c*7r, 
z^O, then 

POO POTT 

Ea,p(z) = / K(a,fi,x,z) dx+/ P(a, p, p, <p, z) d<p, p > 0, p e R, (28) 

Jp J —an 


E a ,p(z) = - 


1 

2 


E a ,p{z) — [ K(a,0,x,z) dx, if P < 1 + a, 

Jo 

sin(aTr) f°° e~*- 1/a 

/ — o 7 — \ . o d A:» J 0=1 + a 

an Jo x l ~ z-X z cos(o:7r) + 2 2 


(29) 

(30) 
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wherein 


K(a t f},x,z) 


x (l-0)/a e - X ^ a 

arc 


X sin[7r(l — /?)] — z sin[7r(l — ft + a)] 
X 2 — 2 x* cos(a7r) + z 2 


P(ot,P, P><P*z) 


pi+(\-P)/a e/ ,1/acos (v , / Q! )[cos(a;) + isin(a;)] 
2onr pe'v -z 


oj = p l ! a sin (p/a) + <p\\ + (1 — P)/oi\. 

The representations in Eqs. (28-30), and similar formulae for the cases | arg z\ = 
an and | arg z | < an presented in Gorenflo et al. [22], are an essential part of 
the numerical algorithm listed below. 


Using the integral representations in Eqs. (26 & 27), it is not too difficult 
to get asymptotic expansions for the Mittag-Leffler function in the complex 
plane. Let a < 2, 0 be an arbitrary number, and (p be choosen to satisfy the 
condition in Eq. (25). Then we have, for any p € N and | 2 r| — > oo: 


(1) whenever |argz| < <p, 


E «A z ) 


z ( e * 1/a 


a 


p z -k 

S r (^ - ak) 


+ 0(|z|- 1-p ). 


(31) 


(2) whenever ip < \ &rgz\ < 7r, 


„-k 


E aA Z ) g T(J3 


ak) 


+ 0{\z\-'->) 


(32) 


These formulae are also used in our numerical algorithm. 


Thorough discussions of properties of the Mittag-Leffler function can be found, 
for example, in Refs. [15,33,45]. 


8.3 Numerical algorithms 


Based on these theoretical results, a numerical scheme listed in Alg. 4 for 
computing the general Mittag-Leffler function E a< p(z ) has been developed 
by Gorenflo et al. [22], which is reproduced here with corrections. Also pre- 
sented is a scheme for computing the derivative of the Mittag-Leffler function 
dE at / 3 (z)/dz, it being Alg. 5, that we use as a memory function in our theory 
of fractional-order viscoelasticity [18]. 
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Algorithm 4 

Computation of the Mittag-Leffler function E a ^{z) + p(z), |/z(z)| < e: 
given a > 0, P € R, z €. C, e>0, 0<C<1 then 
if z = 0 then 

EaA*) = I/W 

elsif a < 1 then 
if \z\ < C then 

k 0 = max{f(l - p)/a 1, “ M)]/ ln (M)l} 

Boj>W = E?=o**/r(/3 + «k) 
elsif \z\ < [10 + 5aJ then 

_ Jmax{l, 2\z\, [- ln(e7r/6)] a }, P > 0 

X0 " [max{(|/?| + 1)«, 2\z\, [-2 ln(*7r/[6(|/?| + 2)(2|/J|)«]r }, 0 < 0 
K(a,0,x,z) = 

u = (f>[l + (1 - P)/a) + p l/a sin ((f) /a) 

P( a,0,P,M = ^P 1+{1 ~ Wa exp[p-/° cos Wa)] 
if | arg z| > ott and 1 1 arg z| - a7r| > e then 
if /3 < 1 + a then 
E a fi{z ) = / O xo K(a, (3, x , «) dX 

else 

^a^(2) = /i X0 #(«, P, X. 2 ) dx + /"*„ P(a, P, 1, 0 , 2 ) d0 

elsif | argz| < o7r and 1 1 argz| - a7r| > e then 
if P < 1 + a then 

E aJ , (z) = /„*« K(a, 0, X ,z) ix+ e- u “ 

else 

EaM = /p|/2 ^(«. P’ x, 2) dx + f-an p ( a , P, \\ z l <£> z ) d 0 + Z z(1 

else 

E a ,p(z) = J£f +1/2 A(a, P, X, z) dx + /“L p («. 0> M + §> 4>, z ) d 0 

else 

fc 0 = L— ln(e)/ln(|z|)J 
if I arg z\ < 3a7r/4 then 

E aJ ,(z) = e* v “ - Efei *-‘/r(P - a*) 

else 

E a A z ) = - Ejtli 2“V r ($ - ak) 

elsif 1 < a < 2 then 

Ea,p{z) = [E a / 2 ,p(\/z) + E a / 2j p(-\fz)\/2 

else 

fc 0 = L«/2J + 1 

E a fi{z) = E^o 1 ^a/fco,^[2 1/fco exp(i27rfc//c 0 )]/A:o 

end 


■«/<* e zl/a 
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Note 2 In Algs. 4 & 5> error tolerance e should be set at e m = machine 
precision, while parameter £ should be set around 0.9 to keep ko reasonably 
sized. 

8. 3. 1 Essentials of the construction 

Algorithm 4 uses the defining series (Eq. 21) for arguments of small magnitude, 
its asymptotic representation (Eqs. 31 &: 32) for arguments of large magnitude, 
and special integral representations (the formulas in Eqs. (28-30) for the case 
where |argz| > air, and similar representations for the cases |argz| = air 
and | arg z\ < an) for intermediate values of the argument that include a 
monotonic part f K(a, /?, x> z) dx an d an oscillatory part f P(a, /?, p, <f>, z) d </>, 
which can themselves be evaluated using standard techniques (cf. App. B). 
These special integral representations are needed because the above series can 
become numerically unstable in this intermediate region where these integrals 
are employed. The outcome is a robust method for computing the Mittag- 
Leffler function. 


8.3.2 Efficiency 

Algorithms 4 & 5 can produce numerical results to any desired level of accu- 
racy. These computations can be expensive, and therefore their use in a finite 
element setting, for example, may be prohibitive. To help meet this need, we 
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Fig. 10. The Mittag-Leffer function E a ^(—x a ) for (3 = 1 and a € (0,2). 


have constructed a table of Pade approximates for E atl (—x a ) in App. C for 
x > 0 and a £ {0.01, 0.02, 0.03, . . . , 0.98, 0.99}. As mentioned above, this form 
of the Mittag-Leffler function arises in many fractional-order models. 


Another scheme based on Eqs. (21 & 32) for solving the Mittag-Leffler function 
E a , i(x) (0.02 < a < 0.98 with a reported relative error that is less than 
1.6 x 10~ 5 ) has been published by Welch et al. [54]. 


8.4 Examples 


Variants of the Mittag-Leffler function E a ^{—x a ) arise naturally in solutions of 
linear, fractional-order, differential equations. In Fig. 10, plots of this function 
are provided for /3 = 1 and a £ (0,2), while in Fig. 11, plots are given for 
a = 1 and 0 £ (0,2). Recall that Ei t i(—x) = e~ x . Parameter a is seen to 
have a strong influence on the overall shape of the curve, while 0 has its most 
pronounced influence on the value of the function at x = 0, in this particular 
case. Numerous example figures can also be found in Gorenflo et al. [22]. 
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Fig. 11. The Mittag-Leffer function E a ^(-x a ) for a = 1 and /} € (0,2). 

9 Summary 


In this paper we provide numerical algorithms that solve fractional-order in- 
tegrals, derivatives of the Caputo type, and differential equations, also of 
the Caputo type. We also present an algorithm that evaluates the Mittag- 
Leffler function and its derivative, which are present in analytic solutions to 
fractional-order differential equations. Collectively, these methods constitute 
a toolbox that engineers can use to address problems where the fractional 
calculus has been used to model some aspect of reality. 
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A Table of Caputo Derivatives 


For the convenience of the reader, we provide this appendix where we give some 
Caputo-type derivatives of certain important functions. 7 We do not strive for 
completeness in any sense, but we do want to give at least the derivatives of 
classic examples. 


Throughout this appendix, a will always denote the order of the Caputo- 
type differential operator under consideration. We shall only consider the case 
ot > 0 and a N, where N := { 1 , 2 , 3, . . .} while N 0 := {0, 1 , 2 ,.. .}. We use the 
ceiling function [a] to denote the smallest integer greater than or equal to a, 
and the floor function [aj to denote the largest integer less than or equal to 
a. Recall that for a € N, the Caputo differential operator coincides with the 
usual differential operator of integer order. 


Moreover, as in Eq. (21), E a t p denotes the Mittag-Leffler function in two pa- 
rameters a and f3 > 0, given by 


E.At) ■= £ 


Sr (ak + py 

while ip is the Digamma function given by 

r'M 


tp(x) := 


rw 


Functions i-Fi and 2 F 1 denote the usual hypergeometric functions, i.e., 


(sometimes called the Kummer, confluent, hypergeometric function [1, Chp. 
13]), the power series being convergent for arbitrary complex z, and 


2 Fi(a,b] c; z) 


F(c) ^ r (q + k)T (b + k) k 

r(a)m^ 0 r(c + k)k\ " 


a, b e R, 


-ci N 0 , 


(the Gauss hypergeometric function [1, Chp. 15]), in which case the power se- 
ries converges for all complex z with \z\ < 1, and may be extended analytically 
into the entire complex plane with a branch cut along the positive real axis 
from +1 to + 00 . (Note that in the formulae below, we only need to evaluate 
this function for negative values of z, so the branch cut for positive z does not 
need to be addressed here.) Finally, i := \f-l is the imaginary unit. 

7 Tables of Riemann-Liouville integrals and derivatives can be found in various 
places in the literature (cf. e.g., Podlubny [45] or Samko et al. [51]). We do not 
repeat those results here. 
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(1) Let f(x ) = a:- 7 . Here we have to distinguish some cases: 


Wf)(x) = 


0 

r U + 1) 


if j € N 0 and j < fa], 

x 3-<* jf j ^ an( j j > J'q-'] or j £ N and j > [aj . 


I T(j + 1 - a) 

(2) Let f(x ) = (x + c) J for arbitrary c > 0 and j € R. Then 

r (j + 1 ) c^'-W- 1 


(/£/)(*) = 


x 


fa ]— q 


T(j-[a\)T(\a]-a + l) 
x 2 Fi( 1, for] - j ; for] - a + 1; -x/c). 


(3) Let /(a:) = x* In a: for some j >L«J. Then 


( d;/)(x) = E(-i) ({) ^£0 Jg j) 


k = 0 


fc J [a] — A: T(j — a + 1) 


+ + 1) xl L«J) -ip(j -<* + l) + lnx). 

(4) Let f{x) = exp(ya:) for some j € R. Then 

(D?f)(x) = 

(5) Let f(x) = sin jx for some j € R. Here again we have two cases: 

( ijWf_1 \\a\/2 \a\-a 

"2f(M - 0 + l) [lFl(1; M " a + 1: «*> 

— i i 7 ! ( 1 ; fa] — a + 1; — ijx)] if fa] is even, 

2 r(M-a + l) l,Fl(1; M - a + 1; ijx) 

+ iFi(l; fa] - a + 1; -ija:)] if fa] is odd. 


Wf){x) = 


(6) Finally we consider / (a:) = cos jx with some j 6 R. As in the previous 
example, we obtain two cases: 


Wf)(x) = 


( j\<A( — \\\fAI2 \a~\-a 

Trad -o + i) llF,(l1 w " a + 1; ijx) 

4- iFi(l; fa] — a + 1; — ija:)] if fa] is even, 

i ? -M('_i)H/2 ;E rai-a 

Tf(R-o + i) M " “ + 1' 

— iFi(l; fa] — a + 1; — ija:)] if fa] is odd. 
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B Automatic Integration 


In Alg. 4 we have need to evaluate two integrals in order to be able to 
compute the Mittag-Lefller function. In this appendix we outline a possi- 
ble strategy for the solution of these integrals. Essentially, we follow the 
ideas explained in quadpack [43]. The routines introduced in that book 
are the generally accepted standard when looking for efficient and reliable 
quadrature algorithms. They are in the public domain, but they can also be 
found in many commercially distributed software packages. The source code 
is written in FORTRAN77 and may be obtained, for example, from the URL 
http://www.netlib.org/quadpack. For integrands of the form encountered 
in Alg. 4, it turns out that the routine DQAG (Double precision Quadrature, 
Adaptive, General purpose) is the method of choice. Our description of this 
routine follows a top-down structure (i.e., we first explain the general strategy 
without giving much information on the details, and at a later stage we fill in 
those details). 


B. 1 Fundamental strategy 


The fundamental idea of automatic integration is the following. The user sup- 
plies the integrand function, the bounds for the interval of integration (i.e., the 
data that define the integral in question), and a desired accuracy (i.e., a bound 
for the relative or absolute error that he/she is willing to accept). The routine 
is then supposed to return either an approximation for the correct value of 
the integral that is sufficiently accurate to satify the user’s requirements, or 
an error flag if it fails to find such an approximation. 

Typically, the algorithms try to achieve this goal by sub-dividing the interval 
of integration in an adaptive way, thus concentrating more quadrature nodes in 
areas where the integrand is difficult to approximate. This leads to a structure 
indicated in Alg. 6. 

In practice, one may also terminate the while loop in this algorithm when, 
for example, too many sub-divisions have taken place thereby using up all 
available memory, or too much computing time has been consumed. In such 
cases, the algorithm should return an error flag. 


B.2 Approximation of the integral 


In Alg. 6 we have left open some of the key details. Specifically, we have not 
said how the approximation of the integral itself will be performed, and we 
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have not indicated how the required error estimates can be found. We now 
turn our attention to these two questions. 


Algorithm 6 

Automatic Integration: 
given a, b G R, / : [a, b] — >• C, e > 0 then 
Calculate an approximation Q[f] for / a 6 f(x ) da: 
and an estimate i for the error / a f(x) da: — Q[f] 

Initialize a list L of the approximations obtained so far 
as L := {([a,b],Q[f],i)} 
while Ei, |e| > e do 

Take the interval from L with the largest error estimate 
and remove it from the list; 

Bisect this interval; 

Calculate approximations and error estimates for the two 
newly obtained subintervals; 

Add these subintervals, their corresponding approximations, 
and their error estimates to the list L 

end 

return the sum of approximations over all the subintervals 


The basic idea behind the solution is the concept of Gauss-Kronrod integra- 
tion. That is, we calculate two different approximations for the integral, both 
of which are of the form 

n 

’yi a j,nf ( x j,n)y 

j = 1 

with suitably chosen values of n, a j>n , and x i<n (j = 1, . . . , n). In particular, we 
begin with a first approximation that is just a Gauss quadrature formula with 
ni nodes, which is the (uniquely determined) quadrature formula that gives 
the exact value of the integral whenever the integrand is a polynomial of degree 
not exceeding 2rq — 1. This formula has been thoroughly investigated. For a 
recent survey we refer to Ref. [3] and the references cited therein. Specifically, 
there is no quadrature formula with the same number of nodes that is exact 
for a larger class of polynomials. Moreover, as stated in [3], both theoretical 
and practical evidence suggest that this method is a very good one. 


B. 3 Approximation of error estimates 


From rather general considerations, it is known that one, single, quadrature 
formula can only give an approximation, but not both an approximation and 
an error estimate. Therefore, to derive also an error estimate, it is necessary to 
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introduce a second quadrature formula. The heuristic argument is as follows. 
Assuming that the second quadrature formula is much more accurate than 
the first, then the difference between the two approximations will be a good 
approximation for the error of the cruder of the two approximations, and 
therefore, a rather conservative (and thus quite reliable) upper bound for 
the error of the finer of the two approximations. Hence, we need a second 
quadrature formula that is significantly better than the first (Gauss-type) 
formula. In view of the quality of the Gauss formula, this can only be achieved 
by selecting a formula with n 2 nodes, where n 2 > n\. In principle, one could 
use a Gauss formula again, but this would be very uneconomical because a 
Gauss formula with ni nodes would have at most one node in common with a 
Gauss formula with n 2 nodes, and so almost none of the information gathered 
so far (i.e., the function values f(xj >ni ), whose calculation is typically the 
most computationally expensive part of the algorithm) could be reused. To 
overcome this difficulty, Kronrod [28,29] suggested to construct a formula that 
is nowadays called the Kronrod extension of the Gauss formula or, shortly, the 
Gauss-Kronrod formula. His formula is based on the Gauss formula with ni 
nodes; it uses n 2 = 2ni + l nodes, and is constructed according to the following 
criteria: 

• The n\ nodes Xj >Tll of the Gauss method form a subset of the n 2 nodes Xj tTl2 
of the Gauss-Kronrod method. 

• The remaining nodes and the weights aj itl2 of the of the Gauss-Kronrod 
scheme are determined in such a way that the resulting method is exact for 
all polynomials of degree not exceeding 3ni + 1. 

A recent survey on Gauss-Kronrod formulas is given in Ref. [13]. 

Algorithms for the concrete calculation of the required nodes and weights 
are available (cf. the survey papers mentioned above and the references cited 
therein). For our purposes, it is sufficient to use the tabulated values given 
in Ref. [43]. In particular, following the suggestions made there, we propose 
using the Gaussian method with 15 points in combination with the 31-point 
Kronrod extension for computing the integral with function K mentioned in 
Alg. 4 (the monotonic part). In view of the nice smoothness properties of 
the integrand, this gives us an approximation with quite high accuracy with- 
out too much computational effort. For the other integral in that algorithm, 
with the oscillatory integrand function P, these formulae may be unable to 
follow the oscillations properly; thus, we propose using the 30-point Gauss 
method together with its 61-point Kronrod extension. This is essentially also 
the method suggested for oscillatory integrals in the NAG Library [40] (see 
the documentation for routine D01AKF). 
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C Table of Pade Approximates for Mittag-Leffler Function 


This appendix presents a scheme for fast computations of 


Ea,i(-x a ), 0 < a < 1, i£l|, 

suitable for use in finite element analyses. As mentioned in §8, this function 
is of prime importance in many applications related to FDE’s. 


The defining series, Eq. (21), is used on the interval 0 < x < 0.1, and its 
asymptotic series, Eq. (32), is used whenever x > 15, as in Alg. 4, but in 
contrast with Alg. 4, a Pade approximate (or rational polynomial) is applied 
inbetween where 0.1 < x < 15; specifically, 


Ea, l(-X a ) 


(~x) ak 

^r(l + ak)' 

a 0 + a,\x + a 2 x 2 
1 + b\X + b 2 x 2 + bzx 3 ’ 
A ( -x)~ Qk 

l £sr(i-o*)’ 


0 < x < 0.1, 
0.1 < x < 15, 
x > 15, 


(C.l) 


where the coefficients a 0} ai, a 2 , bi, b 2 , and b$ are listed in Table C.l for values 
of a varying between 0.01 and 0.99 by increments of 0.01. Also listed are the 
r 2 resultants from each nonlinear regression with the tolerance set sufficient 
to achieve accuracies of r 2 > 0.999, except for a > 0.9 where such tight r 2 
values could not be achieved using the order of approximate employed. The 
Pade approximates were obtained from fits to exact values (within machine 
precision) for the Mittag-Leffler function (obtained from Alg. 4) at 150 evenly 
spaced grid points over the interval [0.1, 15]. The form of this rational polyno- 
mial was selected after a study of numerous tables for function approximation 
listed in the appendices of Hart et al. [25]. Our fits were constrained to be 
exact at the end points so that transitions will be smooth when crossing a 
boundary between solution domains. 

A number written as 0.123456(7) stands for 0.123456 • 10 7 = 1, 234, 560. 


35 


Table C.l: Coefficients for Pade approximates for E a ,i(—x a ), x € [0.1, 15]. 


a 

ao 

ai 

a 2 

b i 

b 2 

63 

r 2 

0.01 

-.424129(3) 

-.695856(6) 

-.391944(6) 

-.138637(7) 

-.795778(6) 

-.207400(3) 

0.999371 

0.02 

-.291078(4) 

-.231705(7) 

-.121393(7) 

-.460100(7) 

-.250385(7) 

-.127923(4) 

0.999345 

0.03 

-.159594(4) 

-.883786(6) 

-.485177(6) 

-.174711(7) 

-.101562(7) 

-.832687(3) 

0.999381 

0.04 

-.401668(4) 

-.151329(7) 

-.678102(6) 

-.298958(7) 

-.144743(7) 

-.138260(4) 

0.999114 

0.05 

-.480619(4) 

-.145550(7) 

-.643881(6) 

-.286638(7) 

-.139751(7) 

-.171118(4) 

0.999129 

0.06 

-.358700(4) 

-.953746(6) 

-.452754(6) 

-.186815(7) 

-.996694(6) 

-.160873(4) 

0.999313 

0.07 

-.431271(4) 

-.958824(6) 

-.424475(6) 

-.187447(7) 

-.952144(6) 

-.175012(4) 

0.999215 

0.08 

-.104813(5) 

-.220562(7) 

-.109649(7) 

-.427837(7) 

-.248747(7) 

-.594850(4) 

0.999401 

0.09 

-.113182(5) 

-.204162(7) 

-.931135(6) 

-.395568(7) 

-.215524(7) 

-.560193(4) 

0.999340 

0.10 

-.100534(5) 

-.165546(7) 

-.753470(6) 

-.319371(7) 

-.177305(7) 

-.527749(4) 

0.999361 

0.11 

-.771890(4) 

-.116952(7) 

-.530686(6) 

-.224658(7) 

-.126954(7) 

-.429168(4) 

0.999382 

0.12 

-.209949(4) 

-.331485(6) 

-.178194(6) 

-.628014(6) 

-.428531(6) 

-.183967(4) 

0.999349 

0.13 

-.201966(5) 

-.264628(7) 

-.118361(7) 

-.503944(7) 

-.292906(7) 

-.123451(5) 

0.999408 

0.14 

-.144185(5) 

-.177380(7) 

-.787251(6) 

-.336281(7) 

-.198184(7) 

-.923408(4) 

0.999420 

0.15 

-.358699(4) 

-.414163(6) 

-.180624(6) 

-.782040(6) 

-.463020(6) 

-.235432(4) 

0.999423 

0.16 

-.367384(5) 

-.387982(7) 

-.157708(7) 

-.732069(7) 

-.414017(7) 

-.218683(5) 

0.999358 

0.17 

-.834348(4) 

-.877610(6) 

-.381391(6) 

-.163991(7) 

-.101120(7) 

-.619007(4) 

0.999453 

0.18 

-.665920(4) 

-.661019(6) 

-.279508(6) 

-.123090(7) 

-.755801(6) 

-.495376(4) 

0.999452 

0.19 

-.880313(4) 

-.839225(6) 

-.352670(6) 

-.155465(7) 

-.970538(6) 

-.689880(4) 

0.999463 

0.20 

-.516937(4) 

-.463982(6) 

-.186660(6) 

-.857478(6) 

-.525412(6) 

-.391939(4) 

0.999446 

0.21 

-.181654(5) 

-.159929(7) 

-.655973(6) 

-.293380(7) 

-.187316(7) 

-.153968(5) 

0.999480 

0.22 

-.287739(5) 

-.241617(7) 

-.962412(6) 

-.441685(7) 

-.280701(7) 

-.244057(5) 

0.999477 

0.23 

-.283969(5) 

-.233523(7) 

-.935683(6) 

-.423970(7) 

-.277398(7) 

-.261404(5) 

0.999497 

0.24 

-.218950(5) 

-.174274(7) 

-.687682(6) 

-.314803(7) 

-.207911(7) 

-.208550(5) 

0.999503 

0.25 

-.184865(5) 

-.142847(7) 

-.555414(6) 

-.256670(7) 

-.171282(7) 

-.182554(5) 

0.999510 

0.26 

-.154930(5) 

-.116572(7) 

-.448241(6) 

-.208289(7) 

-.140906(7) 

-.159998(5) 

0.999520 

0.27 

-.252656(5) 

-.187367(7) 

-.720787(6) 

-.332362(7) 

-.230555(7) 

-.280571(5) 

0.999534 

0.28 

-.348286(5) 

-.250157(7) 

-.938597(6) 

-.441764(7) 

-.306951(7) 

-.392527(5) 

0.999538 

0.29 

-.223178(5) 

-.158819(7) 

-.597404(6) 

-.278265(7) 

-.198723(7) 

-.272201(5) 

0.999553 

0.30 

-.262152(5) 

-.180630(7) 

-.659350(6) 

-.315184(7) 

-.224602(7) 

-.321389(5) 

0.999556 

0.31 

-.110316(5) 

-.741001(6) 

-.263517(6) 

-.128658(7) 

-.919175(6) 

-.137449(5) 

0.999558 

0.32 

-.344887(5) 

-.231540(7) 

-.831618(6) 

-.398374(7) 

-.294522(7) 

-.473795(5) 

0.999578 

0.33 

-.253440(5) 

-.165943(7) 

-.579558(6) 

-.284133(7) 

-.210336(7) 

-.352888(5) 

0.999583 

0.34 

-.320037(5) 

-.206983(7) 

-.712936(6) 

-.351975(7) 

-.264342(7) 

-.467058(5) 

0.999592 

0.35 

-.788648(5) 

-.506620(7) 

-.172863(7) 

-.854865(7) 

-.654223(7) 

-.122026(6) 

0.999604 

0.36 

-.289781(5) 

-.195379(7) 

-.707274(6) 

-.324061(7) 

-.268137(7) 

-.555978(5) 

0.999616 
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Table C.l: Coefficients for Pade approximates for E a ,i(—x a ), x € [0.1, 15]. 


a 

ao 

ai 

a 2 

b i 

£>2 

63 

r 2 

0.37 

-.142051(5) 

-.971770(6) 

-.354821(6) 

-.159385(7) 

-.136569(7) 

-.302386(5) 

0.999617 

0.38 

-.816551(5) 

-.523578(7) 

-.175872(7) 

-.860767(7) 

-.705295(7) 

-.156056(6) 

0.999642 

0.39 

-.575871(5) 

-.348533(7) 

-.107617(7) 

-.574012(7) 

-.450994(7) 

-.991720(5) 

0.999625 

0.40 

-.173808(5) 

-.117277(7) 

-.408516(6) 

-.187857(7) 

-.167833(7) 

-.429558(5) 

0.999651 

0.41 

-.103477(5) 

-.686777(6) 

-.231234(6) 

-.109392(7) 

-.977034(6) 

-.258727(5) 

0.999700 

0.42 

-.195086(5) 

-.129617(7) 

-.429282(6) 

-.204698(7) 

-.185582(7) 

-.514821(5) 

0.999681 

0.43 

-.333165(5) 

-.221252(7) 

-.718492(6) 

-.346552(7) 

-.318194(7) 

-.922518(5) 

0.999693 

0.44 

-.473286(5) 

-.311742(7) 

-.981295(6) 

-.484992(7) 

-.447200(7) 

-.134275(6) 

0.999708 

0.45 

-.147598(5) 

-.982326(6) 

-.304968(6) 

-.151324(7) 

-.142189(7) 

-.447789(5) 

0.999719 

0.46 

-.152914(5) 

-.104712(7) 

-.326140(6) 

-.159246(7) 

-.154553(7) 

-.517302(5) 

0.999725 

0.47 

-.167037(5) 

-.114533(7) 

-.346912(6) 

-.172774(7) 

-.169118(7) 

-.587605(5) 

0.999740 

0.48 

-.142611(5) 

-.991368(6) 

-.293586(6) 

-.148061(7) 

-.147062(7) 

-.531760(5) 

0.999752 

0.49 

-.654706(5) 

-.455618(7) 

-.130849(7) 

-.675094(7) 

-.675481(7) 

-.253173(6) 

0.999766 

0.50 

-.188800(3) 

-.135370(5) 

-.383675(4) 

-.198223(5) 

-.202768(5) 

-.797517(3) 

0.999777 

0.51 

-.516448(5) 

-.373170(7) 

-.102716(7) 

-.541357(7) 

-.559798(7) 

-.228252(6) 

0.999790 

0.52 

-.392154(5) 

-.289687(7) 

-.780114(6) 

-.415703(7) 

-.437101(7) 

-.185908(6) 

0.999801 

0.53 

-.135449(5) 

-.106392(7) 

-.289899(6) 

-.150142(7) 

-.164334(7) 

-.749106(5) 

0.999802 

0.54 

-.953841(4) 

-.748349(6) 

-.193214(6) 

-.104887(7) 

-.114460(7) 

-.530462(5) 

0.999822 

0.55 

-.793386(3) 

-.639333(5) 

-.159875(5) 

-.886439(5) 

-.979799(5) 

-.469973(4) 

0.999835 

0.56 

-.599232(4) 

-.443865(6) 

-.952614(5) 

-.619108(6) 

-.646389(6) 

-.288690(5) 

0.999847 

0.57 

-.337393(5) 

-.286190(7) 

-.663678(6) 

-.388535(7) 

-.438767(7) 

-.223581(6) 

0.999859 

0.58 

-.463960(5) 

-.405168(7) 

-.899289(6) 

-.544254(7) 

-.620515(7) 

-.324201(6) 

0.999871 

0.59 

-.151791(4) 

-.125852(6) 

-.241507(5) 

-.169238(6) 

-.185486(6) 

-.904156(4) 

0.999886 

0.60 

-.502729(4) 

-.464569(6) 

-.919187(5) 

-.611566(6) 

-.705142(6) 

-.377940(5) 

0.999896 

0.61 

-.480103(5) 

-.461103(7) 

-.860143(6) 

-.600297(7) 

-.697998(7) 

-.378472(6) 

0.999907 

0.62 

-.671512(4) 

-.708363(6) 

-.131044(6) 

-.906252(6) 

-.108843(7) 

-.626832(5) 

0.999908 

0.63 

-.312780(5) 

-.328986(7) 

-.533551(6) 

-.418611(7) 

-.494518(7) 

-.268704(6) 

0.999928 

0.64 

-.567954(5) 

-.632239(7) 

-.945394(6) 

-.794766(7) 

-.947107(7) 

-.509518(6) 

0.999937 

0.65 

-.149230(4) 

-.182752(6) 

-.259499(5) 

-.226154(6) 

-.275373(6) 

-.151317(5) 

0.999939 

0.66 

-.433396(5) 

-.549988(7) 

-.664464(6) 

-.674519(7) 

-.815645(7) 

-.407908(6) 

0.999953 

0.67 

-.752321(4) 

-.109682(7) 

-.126030(6) 

-.132094(7) 

-.164296(7) 

-.843703(5) 

0.999948 

0.68 

-.503149(5) 

-.782828(7) 

-.716973(6) 

-.932883(7) 

-.115421(8) 

-.500161(6) 

0.999960 

0.69 

-.113475(5) 

-.204800(7) 

-.158167(6) 

-.240096(7) 

-.301897(7) 

-.117538(6) 

0.999955 

0.70 

-.277604(3) 

-.535397(5) 

-.226639(4) 

-.621910(5) 

-.768192(5) 

-.150446(4) 

0.999976 

0.71 

-.551694(5) 

-.158192(8) 

-.673086(6) 

-.178833(8) 

-.232888(8) 

-.532630(6) 

0.999919 

0.72 

-.310305(5) 

-.102892(8) 

0.663348(5) 

-.115056(8) 

-.146317(8) 

0.202771(6) 

0.99994 
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Table C.l: Coefficients for Pade approximates for E a ,i(—x a ), x G [0.1, 15]. 


a 

ao 

ai 

<X2 

b i 

&2 

b 3 

r 2 

0.73 

-.160766(4) 

-.449254(7) 

-.116371(6) 

-.482429(7) 

-.682061(7) 

-.111373(6) 

0.999719 

0.74 

0.341935(4) 

-.241906(7) 

-.669717(5) 

-.252295(7) 

-.376790(7) 

-.808538(5) 

0.999510 

0.75 

0.173763(4) 

-.103357(7) 

-.513330(6) 

-.106853(7) 

-.212490(7) 

-.787594(6) 

0.999166 

0.76 

-.198004(4) 

-.170209(6) 

-.546410(6) 

-.212224(6) 

-.727938(6) 

-.906762(6) 

0.999483 

0.77 

-.255574(4) 

-.399973(6) 

-.408010(6) 

-.461840(6) 

-.889107(6) 

-.736776(6) 

0.999574 

0.78 

-.107310(4) 

-.102865(6) 

-.104591(6) 

-.125079(6) 

-.207143(6) 

-.205284(6) 

0.999673 

0.79 

-.494450(4) 

-.463767(6) 

-.357502(6) 

-.563528(6) 

-.802633(6) 

-.767703(6) 

0.999707 

0.80 

-.218066(4) 

-.219047(6) 

-.128115(6) 

-.262855(6) 

-.339154(6) 

-.302753(6) 

0.999721 

0.81 

-.534378(4) 

-.539996(6) 

-.258539(6) 

-.645497(6) 

-.766988(6) 

-.675138(6) 

0.999723 

0.82 

-.675467(4) 

-.688012(6) 

-.275764(6) 

-.819312(6) 

-.908437(6) 

-.800006(6) 

0.999715 

0.83 

-.702848(3) 

-.690578(5) 

-.240789(5) 

-.824283(5) 

-.847142(5) 

-.779293(5) 

0.999699 

0.84 

-.609980(4) 

-.583751(6) 

-.176849(6) 

-.697471(6) 

-.667838(6) 

-.642815(6) 

0.999674 

0.85 

-.103053(5) 

-.915547(6) 

-.247199(6) 

-.110278(7) 

-.963517(6) 

-.101397(7) 

0.999639 

0.86 

-.774771(4) 

-.667011(6) 

-.156824(6) 

-.805409(6) 

-.651177(6) 

-.733497(6) 

0.999594 

0.87 

-.838904(4) 

-.693431(6) 

-.141730(6) 

-.840712(6) 

-.624128(6) 

-.763282(6) 

0.999538 

0.88 

-.447750(4) 

-.359964(6) 

-.629983(5) 

-.437447(6) 

-.299094(6) 

-.396021(6) 

0.999468 

0.89 

-.390787(4) 

-.296304(6) 

-.446099(5) 

-.362820(6) 

-.222076(6) 

-.331368(6) 

0.999382 

0.90 

-.724688(4) 

-.509096(6) 

-.654940(5) 

-.630388(6) 

-.335377(6) 

-.584487(6) 

0.999277 

0.91 

-.152262(5) 

-.100158(7) 

-.106340(6) 

-.125254(7) 

-.574996(6) 

-.117309(7) 

0.999150 

0.92 

-.150108(5) 

-.942028(6) 

-.792810(5) 

-.118669(7) 

-.467003(6) 

-.112366(7) 

0.998998 

0.93 

-.106495(5) 

-.626399(6) 

-.398859(5) 

-.797901(6) 

-.255226(6) 

-766267(6) 

0.998816 

0.94 

-.153906(5) 

-.847297(6) 

-371308(5) 

-.109242(7) 

-.268144(6) 

-.106481(7) 

0.998600 

0.95 

-.600332(4) 

-.311566(6) 

0.756057(4) 

-.406277(6) 

-.710353(5) 

-.401715(6) 

0.998343 

0.96 

-.445797(5) 

-.214492(7) 

-.119855(5) 

-.284119(7) 

-.270455(6) 

-.285641(7) 

0.998041 

0.97 

-.355171(5) 

-.159674(7) 

0.206708(5) 

-.214673(7) 

-.423690(5) 

-.219239(7) 

0.997685 

0.98 

-.145397(5) 

-.610793(6) 

0.190402(5) 

-.834188(6) 

0.464083(5) 

-.865369(6) 

0.997270 

0.99 

-.105996(5) 

-.416113(6) 

0.204340(5) 

-.577797(6) 

0.755603(5) 

-.608788(6) 

0.996786 
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